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We develop a Ginzburg-Landau theory for the Jaynes-Cummings-Hubbard model which effec- 
tively describes both static and dynamic properties of photons evolving in a cubic lattice of cavities, 
each filled with a two-level atom. To this end we calculate the effective action to first-order in 
the hopping parameter. Within a Landau description of a spatially and temporally constant or- 
der parameter we calculate the finite-temperature mean-field quantum phase boundary between a 
Mott insulating and a superfluid phase of polaritons. Furthermore, within the Ginzburg-Landau 
description of a spatio-temporal varying order parameter we determine the excitation spectra in 
both phases and, in particular, the sound velocity of light in the superfluid phase. 



I. INTRODUCTION 

In many fields of physics, especially in the fields of in- 
formation processing, material science, superfluidity, and 
the relatively new field of quantum information theory, a 
profound understanding of strongly correlated quantum 
many-body systems is of striking importance in order to 
further improve existing applications and invent new ones 
[TJ [2] . This is due to the fact that these research fields 
mainly use solid-state systems in which strongly corre- 
lated systems appear quite naturally. However, it is ex- 
perimentally challenging to access the microscopic prop- 
erties of such systems due to the short time and length 
scales involved. Therefore, motivated by Feynman's con- 
jecture of the quantum simulator [3J, artificial structures 
have been considered to create effective many-body sys- 
tems, which can be investigated under much more con- 
trollable and tunable experimental conditions. 



A. Optical Lattices 

The first attempts to build up artificial many-body 
structures used Josephson junction arrays [H [5], which 
proved to be capable of simulating the dynamic proper- 
ties of Bose-Hubbard systems including the generic quan- 
tum phase transition of this model [5] . Additionally, over 
the last two decades, the advances in preparing and con- 
trolling ultra-cold atoms led to new experimental real- 
izations, which raised a huge amount of interest and re- 
search in this field. For example one investigated the 
interference of BEC clouds [3 [8], studied rotating BEC's 
[SHU], observed spinor condensates [T3J, where BEC oc- 
curs in different hyperfine states, analyzed Bose- Fermi 
mixtures [141 I15j . where a pure BEC is contaminated 
with fermions, described BEC's in disordered potentials 
[IMS], dipolar BEC's QH120], and, more recently, tried 
to probe the properties of BEC's in zero gravity |21) . 
Of particular interest for the simulation of strongly cor- 
related quantum many-body systems has been the re- 
alization of BEC's trapped in optical lattices [351 H3]- 



Since their experimental realization, optical-lattice sys- 
tems have initiated intensive studies and led to a multi- 
tude of new applications such as entanglement of atoms 
EH HI], quantum teleportation [25] , Bell state exper- 
iments [27], disorder [28H3T] . and ultra-cold molecules 
[32l [33] . to name but a few. Unfortunately, the experi- 
mental approaches discussed so far face some crucial limi- 
tations. On the one hand, it is necessary to cool down the 
considered system to some nano Kelvin above absolute 
zero and, on the other hand, it is experimentally challeng- 
ing to control and access single sites individually, How- 
ever, recently developed experimental techniques also al- 
low for single site addressability [34H37] . Nevertheless, 
these experiments still need ultra-cold temperatures, a 
restriction which could be circumvented by using cavity 
lattices. 



B. Cavity QED Lattices 

Encouraged by the latest progress in the fabrication 
and manipulation of micro cavities [38 - 40J , Philippe 
Grangier and others [4TH4"f5] proposed a new experimen- 
tal setup using cavity quantum electrodynamics (QED) 
schemes. The underlying idea behind this new approach 
is to build up a lattice from micro cavities and place 
some real or artificial atoms in each cavity, for exam- 
ple Josephson junctions or quantum dots. Subsequently, 
light is coupled into the system in such a way, that it 
interacts with the atoms. As a result, the coupling be- 
tween the light field and the atoms leads to the forma- 
tion of quasi-particles, so called polaritons. These quasi- 
particles behave like real bosonic particles on the lat- 
tice. In fact, Bose-Einstein condensation of polaritons 
was recently experimentally achieved in semi-conductor 
cavities filled with quantum wells |47H49j and even su- 
perfluidity could be observed [50] [51]. This new idea 
for a quantum simulator based on cavity QED does not 
share the above mentioned limitations of the optical lat- 
tice approach. Due to relatively huge distances between 
the cavities, local control and accessibility emerges quite 
naturally for these systems. Hence, it is possible to ana- 
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lyze these systems without destroying them, in contrast 
to the time-of-fhght imaging technique used for ultracold 
quantum gases. Since the atoms are trapped inside the 
cavities right from the start and their thermal motion 
does not quantitatively disturb the polariton dynamics 
[52] , BEC experiments with cavity QED setups can thus 
be performed even at room temperatures |53| . However, 
in order to facilitate stable experiments with this setup, 
one needs a strong coupling between light and matter in 
order to reduce the losses induced by the spontaneous 
emission. Fortunately, over the past few years, this so 
called strong-coupling regime [541 1551 162j has become ex- 
perimentally accessible for a large number of different 
setups [531 (56H6U [63] . 

However, a drawback of the strong coupling is the short 
polariton lifetime which prevents to reach thermal equi- 
librium. Thus, in a real experiment the polaritons decay 
faster then they can equilibrate via phonon emission and 
polariton-polariton scattering and therefore the system 
needs to be continuously pumped by an external laser. 
Recently, F. Nissen et al. theoretically investigated the 
coherently driven and dissipative JCHM [64] . They could 
show that the important photon-blockade effect prevails 
in this system for weak hopping. This effect guarantees 
a strong polariton-polariton interaction which allows the 
system to reach a quasi thermal equilibrium [UJ [3H] . In 
principal one can modify the method presented in this 
paper to include external baths in order to account for 
this non-equilibrium situation. Subsequently tracing out 
the bath degrees of freedom yields an effective theory for 
the reduced system. Nevertheless, we will in this paper 
focus on the quasi equilibrium situation with a steady 
state polariton density in order to demonstrate the the- 
oretical feasibility. 



C. Outline of the Paper 

This paper is structured as follows. In Section 2 we 
will introduce the underlying Hamiltonian for the Jaynes- 
Cummings Hubbard (JCH) model and briefly analyze its 
dynamical properties in some special limits. In Section 3 
we will derive an effective Ginzburg-Landau action for 
the model. To this end we will use the approach of 
Refs. [65, 66J, that has already been successfully applied 
to analyze collective excitations of Bose gases in optical 
lattices [H7[ EE], and transfer it to a cavity QED sys- 
tem. In contrast to other methods, this procedure yields 
a consistent thermodynamic theory for hnite tempera- 
tures in the whole quantum phase diagram. Therefore 
our Ginzburg-Landau approach is of relevance in order 
to develop a thermometer for the JCHM. For instance, 
to establish a thermometer for the usual Bose-Hubbard 
model turned out to be a difficult task [5S] ED] . The re- 
spective results for both static and dynamic properties 
as well as their temperature dependence obtained from 
the effective action will be presented in Section 4. In 
particular we focus on the excitation spectra and on the 



sound velocity of light and discuss their dependence on 
the experimentally accessible detuning parameter. 

II. THE MODEL 

For our lattice model we consider a Hamiltonian of the 
form 

H = H Q +H h . (1) 

Using the convention h = 1 throughout this paper, the 
local part of this Hamiltonian is given by 

Ho = [( w - ^fii + A afar + g (d, d+ + at drj 

i 

(2) 

and consists of a sum over local Jaynes-Cummings Hamil- 
tonians in the rotating wave approximation (RWA) [71j . 
Here u denotes the frequency of the local monochromatic 
photon fields associated with the annihilation and cre- 
ation operators dj and dt. The parameter A = ui — e 
corresponds to the detuning between the local mode fre- 
quency ui and the energy splitting e of each two-level 
system. The operators and a~ are the ladder oper- 
ators of the two-level systems and g is the real coupling 
strength between each two-level system and the bosonic 
mode. The operator hi = djdj + d^d^ represents the 
on-site polariton number operator and N = J^- hi is the 
total polariton number operator. As it commutes with 
H n , the number of these polaritons, which are coupled ex- 
citations of the two-level system and the local light fields, 
is a conserved quantity in this model. Finally, working in 
the grand-canonical ensemble, yields an additional term 
in the Hamiltonian which is proportional to the polariton 
number operator and the chemical potential \i. 
The Jaynes-Cummings model is well known in the lit- 
erature and the underlying Hamiltonian can be exactly 
diagonalized within the RWA [75] leading to the energy 
eigenvalues 

E n± = -n eS n+ i [A±R n (A)} , n > 1, (3) 
E = 0, n = 0. (4) 

Here, we introduced the generalized Rabi frequency 
R n (A) = ^A 2 + 4 g 2 n and the effective chemical po- 
tential /i ff = /i — ui. The eigenvalues E n ± correspond to 
the energy eigenstates 

\n, +) = sin#„ |n, g) + cos(9„ \n — 1, e) , (5a) 
\n, —) = cos#„ \n,g) — sm8 n \n — 1, e) , (5b) 

with the mixing angle 9 n = | arctan ^ 2g ^" ^ . The vac- 
uum state is given by |0) = |0,<?). We note that the 
energy spectrum for a fixed polariton number n > nat- 
urally splits into an upper and a lower branch, where the 
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lower branch is always lower in energy than the upper 
branch. 

The possibility for photons to tunnel between next neigh- 
boring cavities is modeled by a Hubbard-like hopping 
Hamiltonian of the form 



This Hamiltonian can be diagonalized in Fourier space 
leading to 



(9) 



Hi 



I, — ^ 

<i,3> 



oja,- 



where the sum runs over all next neighbor lattice sites. 
As the tunnel matrix elements Kij exponentially decay 
with increasing distance between the lattice sites i and j, 
it is justified to assume that if i and j are next 

neighbors and Kij = otherwise. 

As a first analysis of the ground-state of the model de- 
scribed by Eq. (JlJ we consider the physically relevant 
extremes of both the atomic and the hopping limit. In 
the atomic limit k -C g, Eq. ([I]) simplifies to H w Hq 
which decomposes into purely local contributions with 
eigenvalues ([3]). Obviously, in this regime the photons 
can not move in the lattice and, thus, all excitations are 
pinned to their respective lattice sites. Therefore, the 
ground state wave function of the whole lattice is simply 
a direct product of the local on-site ground state wave 
functions. For this reason the ground state of the whole 
system is reached when each on-site Jaynes-Cummings 
system is in the lowest energy state Eq or one of the 
lower-branch states E n _. However, decreasing the differ- 
ence uj — fi one eventually reaches a point when Eq = -El- 
and, hence, adding a polariton excitation becomes ener- 
getically favorable. The successive repetition of this ar- 
gument leads to a complete set of such degeneracy points 
E n — = Er n+ i^_, which are characterized by the explicit 
relations: 



^ = -L [R n (A) - R n+1 ( A)] , n > 1 (7a) 
9 2# 



(6) with the energy dispersion e(k) = — /x e ff — 2 k ^2 cos(fc i a), 

i=l 

where a is the lattice constant of a simple d-dimensional 
cubic lattice. Thus, the Hamiltonian of the hopping limit 
^ is local in Fourier space. This situation corresponds 
to the superfluid phase of the system. Signatures of this 
phase have already been observed in the interference pat- 
terns of time-of-flight experiments with Bose-Hubbard 
systems |74) . 



= J- [A -IMA)] 

9 2g 



n = 0. (7b) 



From the above discussion follows that, in the atomic 
limit, the local polariton number is fixed at each lattice 
site for a given set of parameters. When the total polari- 
ton number equals an integer multiple of the number of 
cavities, this regime is called the Mott insulating phase 
which has been predicted [6] and experimentally observed 
in the Bose-Hubbard model [73]. 

As a second interesting limit we investigate the ground- 
state in the regime when k 3> g, i.e. where the photon 
hopping dominates the system dynamics. We addition- 
ally assume that all two-level systems are in their respec- 
tive ground-state as the system minimizes the energy. 
Hence, we can drop all atomic contributions, which leads 
to the hopping-limit Hamiltonian 



H 



"Meff 



N ■ 



<i,j> 



(8) 



III. EFFECTIVE ACTION 

In the following section we derive an effective action for 
the JCH model from the free energy. Explicitly calculat- 
ing the lowest hopping order contributions of the effective 
action amounts effectively to a resummation of infinitely 
many hopping contributions. Therefore, this effective ac- 
tion allows us to determine the quantum phase transition 
as well as calculate the excitation spectra, energy gap, ef- 
fective mass, and sound velocity for finite temperatures 
in the Mott phase and in the superfluid phase. 



A. Free Energy 

Following the approach used for example in Refs. 
l66| we additionally introduce source currents j(r), j*(r), 
in the system Hamiltonian ([I]) leading to the new Hamil- 
tonian 



H'(r) [j, j*] =H + Y, [fi ( T ) ^ + 3i (r) a\ 



(10) 



which is now a functional of the currents. These artificial 
currents, which explicitly depend on the imaginary-time 
variable r, will be used to artificially break the U{1) sym- 
metry that is responsible for the quantum phase transi- 
tion in the model [75]. Since all physical results are ob- 
tained in the limit of vanishing currents, we can treat all 
terms in the Hamiltonian (10), which are proportional 
to the source currents, as small quantities. Furthermore, 
for small hopping amplitudes Kij all off-diagonal contri- 
butions in Eq. ( 10 ) become small. Hence, we decompose 



the new system Hamiltonian into the form 
H'(r)\j,r]=H +H 1 (r)\j,f] 



(11) 



where the local part Ho from Eq. ^ is exactly solvable 
and the remaining part 



#i(r) [j,f] =i?h + E [#( r ) & i + S l 



(12) 
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with the hopping Hamiltonian from Eq. ^ can be 
treated as a perturbation. 

We aim at establishing a thermodynamic perturbation 
theory in the present section. Therefore, it is convenient 
to switch from the Schrodinger picture to the imaginary- 
time Dirac interaction picture. This leads to a refor- 
mulation of the partition function as a perturbation se- 



ries, involving just the perturbative part ( 12 ) of the full 
Hamiltonian of the system. Introducing the abbrevia- 
tion (») = -g^Tr|« e _/3 ^°| with the inverse tempera- 
ture /3 = l/(fceT) for the thermal average with respect 
to the unperturbed system, the partition function takes 
on the following form 



Z = Z (U D ((3,0) 



Here, the partition function of the unperturbed system 
is given by 



Z = Tr 



(14) 



and the imaginary-time evolution operator in the Dirac 
picture is defined as 



thus have to calculate the power series up to fourth order 
in Hence, we expand the logarithm in expression 

( |17[ ) and keep all terms up to fourth order in j and j* 
and first order in k. 

This procedure leads to an expansion of the free energy 
functional in terms of imaginary time integrals over sums 
of products of thermal Green functions with respect to 
the unperturbed system of increasing order. The nth 
order thermal Green function with respect to the unper- 
turbed system is defined as 



( r lJ*li 



T 



4 K) < (n) ...a], «) a tn (r n ) 



(18) 



(13) In principle, one could now make use of the definition 



(18) and calculate the expansion coefficients of the free 



energy straightforwardly. However, with increasing order 
of the thermal Green function the calculation becomes 
more and more complex due to the increasing number of 
space- and time-index permutations. Therefore, we use 
another approach to calculate the thermal Green func- 
tions, which automatically takes care of the emerging 
problems. 



C/bGM) = ^exp I -J dTH x {r) \j,f] V . (15) 



B. Cumulant Expansion 



Note that all operators, which depend on imaginary-time 
variables, have to be taken in the imaginary-time Dirac 
interaction picture, i.e. Ot>(t) 



,H t Qg—H a t 



Using the above definition (12 1 together with Eqs. (13) 



and (151 we see that the partition function also becomes 
a functional of j(r) and Splitting the grand- 

canonical partition functional into the respective pertur- 
bative contributions 



Z\j,j*]=Z \l + J2Z n \j,; 



(16) 



the free energy defined by T = —(3 1 In Z can be written 
as 



(17) 



71=1 



Here we introduced the free energy of the unperturbed 
system as the usual expression Fq = — (3~ 1 \nZ . In 
the next step we expand the free energy functional in 
a power series of the perturbation parameters j,j*, Ku. 
In this paper we will focus just on the lowest order con- 
tributions from the hopping and thus neglect all terms 
of higher than first order in Kij. Furthermore, accord- 
ing to the Landau theory [73], one needs to consider all 
terms at least up to fourth order in the order parameter 
to describe the thermodynamic properties of a second or- 
der phase transition. Since we will see later on that the 
source currents are of the order of the Ginzburg- 
Landau order parameter for the considered system, we 



Usually one would apply in field theory the Wick theo- 
rem to decompose n-point correlation functions into sums 
of products of 2-point correlation functions [77] . Unfor- 
tunately, this is not possible for the considered system, 
since the Wick theorem just holds for systems, where 
the unperturbed Hamiltonian is linear in the occupation 
number operator. Instead, in our case one has to use 
the so called cumulant expansion, which was originally 
developed for the Hubbard model as is reviewed by Met- 
zner [75]. It states that the logarithm of the partition 
function is given by the sum of all connected Green func- 
tions. The power of this approach lies in the fact that 
these connected Green functions can subsequently be de- 
rived from a single generating functional by performing 
functional derivatives with respect to the currents. Due 
to the fact that the unperturbed Hamiltonian ([2| decom- 
poses into a sum over local contributions, the generating 
functional decomposes into products of purely local cu- 
mulants: 



-fdr[j^T)al(r)+j*(r)aAr)y 



The local cumulants then follow from 



cl 0) («; 



Tn) 



(19) 



(20) 



5j i>i (r[)...5j Vn (r n )Sjt 1 (n)...5jtjT n ) 



j=j*=0 
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Due to the local structure of ( 19 



ish unless the site indexes are a 
relation 



all cumulants ( 20 1 van- 



1 equal which yields the 



(*1 ) r l i • • ■ i *n i T n 1*1 ) T l ! ■ • ■ i *1 i T n 



(21) 



Hence, we just have to calculate the local cumulants 
Cn (ii\T~[, . . . ,t„\ti, . . . ,T n ). Performing the calcula- 



tions according to formula ( 20 ) and rearranging the re- 



sulting terms yields the cumulant decomposition for each 



thermal Green function (18). The lowest order Green 
functions read 



(i, n \j, t 2 ) = c[ 0) (i;n|r 2 )^, 

(0) 



G 2 0) (hn;j,T2\k,T 3 ;l,T4) = (i;Ti,T 2 |T 3 ,T 4 ) SijSjkS, 



+C[ 0) (i;n\T 3 ) CT (i;r 2 |T 4 )^ fe ^ 



(0) 



+C{ 0) (i;ri|r 4 ) C x w (j;r 2 \r 3 ) 6 i{ 6 jk 



(o) 



(22) 



With this cumulant decomposition we find the following 
expansion of the free energy functional 



1 f p f p \ 

F [j, f] = Fo ~ ~b J dTl J dl ~ 2 { 

1;3 > 



4 0) (i;n\T2)Sij + a { 2 ) (i;n\T 2 ;j) 



i ft 

+ t / dr 3 dr 4 
4 Jo Jo 



a\ («';ti,t 3 |t 2 ,t 4 ) S i3 + 2a\ (i;Ti,T 3 |r 2 ,T 4 ;j) 



ii(Ti)i*(r 2 ) 

ii(n)ii(T 3 ) j*(T 2 )j*(r 4 ) 



i /"' 3 /" m I 

+ 2/ dr 3 / dr 4 fi 4 (*;n,T 3 |r 2 ,r 4 ; j) jj(n) ji(r 3 ) jffo) jj(r 4 ) L 



(23) 



where the introduced expansion coefficients are defined 
as 



4 0) (*;n|r 2 ) = cf (i; n |7a) 



(24) 



a 2 ±J (*;n|T 2 ;j) = «ij / drC 1 0) (i;n|r) Ci u; (j;t\t 2 ) , 
Jo 

<4 (*; t i,t- 3 |t 2 ,t 4 ) = cf ] (i;n,T 3 \T2,T 4 ) , 

rP 

{i;T 1 ,T 2 \T 3 ,T 4 ;j) = Kij / dr Cf ] (i; r, r 2 |r 3 , t 4 ) 
Jo 

x^i 0) (j;n|r), 

4 1 (*;n,T 3 |r 2 ,r 4 ; j) = Kij / dTC^ 0) (i;ri,T 2 |r 3 ,r) 

Jo 

x CP (j;r|r 4 ). 

Due to the locality of H in Eq. ([2| the cumulants ci°^ 
do not depend on the site indexes For this reason, we 
will drop the site index in the following calculations for 
convenience. Finally, we notice that the form of the above 
coefficients ,a 4 , and a 4 can be further simplified 
by going into frequency space. Therefore we perform the 
Matsubara transformation 



fM = ^=^drf(r)e l 

1 OO 



(25) 
(26) 



with the Matsubara frequencies 
2 7r m 







m £ Z. 



(27) 



At first, we calculate the coefficient a 2 °^ (w m i |w ro2 ) in 
Matsubara space. Due to frequency conservation the fol- 
lowing relation has to hold 



<4 (wmi|w m2 ) 



a 2 0) Ki) S u 



(28) 



This coefficient can be derived from the expression (24) 



with the help of relations (18) and (22) by performing 
a Matsubara transformation 
mapping introduced in Ref. 



251) . Using the polariton 
to calculate the thermal 



expectation values, we obtain the following result: 



,(0) 




71=1 



-PE no 



(*(n+l 



E na — S(n+l)a' + * w n 



(n-l)c 



-En 



(29) 



The coefficients t na p in the above expression stem from 
the fact that there exist two kind of polariton species, 
where the lower branch is labeled by a,/3 = —1 and the 
upper branch by a, /3 = +1. These coefficients are de- 
fined as 



ra— — 00 



t n ±- 



fa a n K- 1 + y/n-l b± 6„_ : , (30a) 



at-! + Vn - 1 b* a n _ 4 , 



(30b) 
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with mixing angle dependent amplitudes given by 



sin6» n , a = + a 
" s cos0 n , a = - 



cos 6 n , a = 
— sin 9 n , a 



(31) 

With the help of this result we can also determine the 
higher hopping corrections. By using frequency conser- 
vation again we find the relation 

(wmi |w m a) =a 2 0) (w m i) a 2 0) (w m2 ) S Uml>Uma (32) 



and 



(0)/ I \ _ 1 r 

{- 4 0) ( w mi) af ] (u) m3 ) [S UmuUm3 S Um3tUmi 
i f 



x (f[tf(n)a<(T 3 )a(T2)a{T i )] 



x e 



i ( — Ii) a l Ti + w m2 T 2 — W m 3 T3 +W m4 T 4 ) 



} 



(33) 



The calculation of the latter expression is complicated 
and rather lengthy. Therefore, we put a detailed cal- 
culation of this quantity in the appendix. Nevertheless, 
from general considerations like frequency conservation 
and integral properties in Matsubara space, we can de- 
duce right away that the first order hopping correction is 
of the form 



(!) i i \ _ C°) i \ 

«4 ^ml,<^m3|W m 2,U; m 4j — &2 l w m2 j 
X 4 (^ml 5 W m3 |w m 4)<W 



(34) 



C. Ginzburg-Landau Theory 



Within this section we finally derive the Ginzburg- 
Landau action for the Jaynes-Cummings-Hubbard 
model, which is the proper thermodynamic potential to 
describe the quantum phase transition of this system. 
Since the symmetry-breaking currents are no phys- 
ical quantities, one has to transform them into physical 
fields in terms of the order parameter. This is accom- 
plished by means of a Legendre transformation of the 
free-energy functional to an effective action as has al- 
ready been explored in detail in the context of the ther- 
mal phase transitions [80j HI]. In order to do this in 
a concise way, we first rewrite the grand-canonical free 
energy ( 23 ) in the following form in Matsubara space 



F [i.J*] = -Ft) - o Yl Yj ( W ml, W m2 ) (Uml) 

i,j W m i,W m 2 
k,l WHi3,W m 4 

X ji(Wml)ji(w m 3)ifc(w J n2)j*(Wm4)} + ■ ■ ■ , (35) 



where we have introduced the abbreviations 
Mij (w m i ,u m2 ) = a 2 0) (wmi ) 



+ Kij (w m l) «2 '' ( W m2) 



(0) 



J W m l,W m 2 5 



(36) 



and 



Nijkl (w m i, W m 2, W m 3, W m 4) — -a 4 ° (w m l, W m 3|w m 4) 
^ij" #fc,i <^i,fc + K ik (Wma) <5j,; 



+ (w m3 ) 4,; 



-'lil m i +!il m3 ,tu m 2 t"i 



(37) 



Now we define the Legendre transformation by self- 
consistently introducing the order parameter field 
^(wm) according to 

%(u, m ) = (a i (LJ m )) =i3 (38) 

Note, that this Ginzburg-Landau order parameter field 
differs from the Landau order parameter by being space 
and time dependent. Inserting expression (35 1 into equa- 



tion (38) yields the following relation for the Ginzburg- 



Landau order parameter field 



*i (w m ) = ~ y^^~] S M pi (w m i,w m ) jp (w ml ) 

— ^Ipki (w m i, W m 2, W m 3, W m ) 

x ji(w m i)jp(w m3 )jfc(a; m 2)} + . . . . (39) 



Furthermore, relation |38|) motivates to introduce the ef- 
fective action 



/3 kiu<]ii 



(40) 



where VP and j are conjugate variables satisfying the Leg- 
endre relations 

3*M = -b J T ^ . , j*(0=-/^ ^— ■ (41) 



<5* 4 K) 



Using the fact, that physical situations correspond to 
vanishing currents j = j* = 0, yields the following equa- 
tions of motion 



ST 



ST 



SV* (wj S^i (w m ) 



= 0. 



(42) 



Thus, the effective action is stationary with respect to 
the order parameter field. Now, in order to determine the 
explicit form of the effective action, we need to express all 
symmetry-breaking currents j by the Ginzburg-Landau 
order parameter field Therefore, we recursively invert 
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relation (39) up to first order in the hopping strength n, 



which yields 

ji (U m ) = ~ E E M iV (^m, Wml) 

— 2 E E Ni q kp (w m l , W m 2 , W m 3 j W m ) 

CJ m 2 ,o; m 3 

X Jl(u> m l)J q (u m3 ) Jfc (^ntf)] + ■ ■ • , (43) 

where we define the abbreviations 

■h W = - E E ' Wm ) *p ) , (44) 



P "ml 



and 



M^ 1 (w m l,U> m2 ) 



4° (i.^ml) 



(°) C • 1 



(45) 

Inserting the above relations into the Legendre transfor- 
mation ( 40 1 gives the explicit expression for the effective 



action up to the desired order in the hopping strength 

2 



«4 ) (^ml,^m3|^m2,^m4) jPj (^mg) g| (<^m4) 1 

2 0) (wmi) a 2 0) (w m2 ) a^ 0) (w m3 ) a^ 0) (w m4 ) J 



4 a 



(46) 



Equation (46) is the thermodynamic potential for the 



Jaynes-Cummings Hubbard model up to the desired ac- 
curacy in both the order parameter and the hopping pa- 
rameter. However, as should be clear from the approx- 
imations performed so far, this expression can, in prin- 
ciple, be extended to include higher order corrections. 
In the next section this result is used to analyze several 
properties of the considered system. 



IV. RESULTS 

Having derived the effective action in the previous 
chapter, we now use this result to extract both thermody- 
namic and dynamic properties of the Jaynes-Cummings- 
Hubbard model. The starting point for this analysis is 
the equation of motion (42). Inserting (46) yields 



= 



1 



4 (^m) 



7 , K ij 



E 

W m l ,UJ m 2 ,^m3 



(47) 



x <4 0) (q; m i,^ m3 |u; m2 ,a; m ) ^ (a; m i) ^ K3) jg| K2) 

2 a 2 0) (uj m i) 4 (w m2 ) (w m3 ) a^ 0) (w m ) 

Using a particular ansatz for the order parameter field 
allows to examine both static and dynamic order param- 
eter fields (uj m ). 



A. Static Results 

First, we consider an equilibrium situation, where the 
order parameter field is constant in both space and time 

*i(w m )=v / l8*r<J«--,o- (48) 



Inserting this ansatz in the stationarity condition (47) 
yields the following relation for the equilibrium order pa- 
rameter 



,(°) 



(0) 



P 4 0) (0,0|0,0) 



,(0) 



(0) KZ 



(49) 



where z = 2d represents the coordination number of the 
(i-dimensional cubic lattice. Since the order parameter 
field is zero in the Mott insulator phase and takes on 
finite values in the superfluid regime, we can extract the 
quantum phase boundary from the condition that the 
equilibrium order parameter ( 49 ) has to vanish. Thus 



the quantum phase boundary is defined via the relation 



1 



(0) 



(50) 



Together with the result ( 29 ) this equation yields a phase 
diagram which is pictured in Fig. [T]for vanishing detun- 
ing A = 0. Here we plot the effective hopping strength 
versus the effective chemical potential leading to a lobe 
structure where each lobe is associated to a specific mean 
on-site polariton number. The regions within these lobes 
correspond to the Mott insulator phase, whereas the ex- 
terior region corresponds to the superfluid regime. First 
we note that the phase boundary for zero temperature 
is consistent with the results from Refs. [751 IH2]- More- 
over, our Ginzburg-Landau theory also yields the phase 
boundary for finite temperatures. From Fig. [I] we can see 
that increasing the temperature leads to a smeared out 
phase boundary. Thermal fluctuations mostly affect the 
region between two neighboring Mott lobes, whereas the 
middle of each Mott lobe almost does not change. This 
effect is stronger for lobes with higher polariton number 
as these system configurations are rather unstable. Since 
a genuine Mott insulator is defined by a vanishing com- 



pressibility Kt 



1 



dn 2 



we find that the Mott 



insulator phase is only present for zero temperature. For 
finite temperature it becomes a mixture of normal and 
Mott insulating phase instead [53]. 

Note that, in order to validate our results obtained so 
far, we follow Ref. |65| and explicitly compare our ex- 



pressions for a 2 
calculation. We 

J~MF — J~0 — J2i 

identifying ji(r) 



(0) JO) 



get 

„MF 



with a corresponding mean-field 
back the mean-field free energy 
|* 4 | 2 + i/3ari^| 4 ) by formally 



and 



_MF 



4°'(0,0|0,0) (kz) 4 . With this we reproduce 
the mean-field results of Refs. [79l |82] including the de- 
pendence of the quantum phase boundary on the detun- 
ing parameter A. If the system is tuned out of resonance, 



-kz^u af F = 4 0) (0)(kz) 2 



K Z 




FIG. 1. (Color online) Quantum phase boundary for zero 
and finite temperatures at resonance A = 0. Interior of the 
lobes for T = OK corresponds to the Mott insulator phase 
whereas the exterior corresponds to the superfiuid phase. For 
finite temperatures the Mott insulator vanishes and becomes 
a mixture of normal and Mott phase. Note that the drop 
of the phase boundary at fi e g ~ —O.lg is a numerical rem- 
nant resulting from a summation cutoff in Eq. (29 1. A full 



summation would lead to an infinite sequence of Mott lobes 
approaching jU e a = 0. 



i.e. A / 0, all Mott lobes with higher mean particle 
number than one shrink in size and are shifted to smaller 
values of the effective chemical potential [x e g irrespective 
of the sign of the detuning. The Mott-lobe with mean 
particle number 1 grows in size for negative detuning 
and shrinks for positive detuning, whereas the Mott-lobe 
with mean particle number covers the rest of the phase 
diagram and, thus, is the only Mott insulator region that 
does not form a closed lobe. This special behavior of the 
first two Mott regions stems from the composite nature of 
the polaritons and has no analogue in the Bose-Hubbard 
model. 



B. Dynamic Results 

Within this section we analyze the dynamic behavior 
of the JCH model. We especially focus on signatures 
within excitation spectra such as energy gaps, effective 
masses and the sound velocity of polaritons. These sys- 
tem properties are experimentally accessible via emis- 
sion and transmission spectroscopy |84j .In order to derive 
these properties we investigate the dynamic behavior of 
the effective action around the equilibrium fields. Thus, 
introducing the vector = (ty, $*) we insert the ansatz 
* = * 



eq 



in the equation of motion (47) and its 



conjugate complex. This yields a system of two coupled 
equations corresponding to the equations of motion for 
the elongations around the order field. In order to 
find a non-trivial solution of the equation of motion, we 



obtain the relation 

S 2 T S 2 T 

<s* 4 ( Wm ) (-w m ) st>* (-uj m ) m*j ko 
<5 2 r s 2 r 



o = 



(51) 



<S*i (w m ) SV* (u> m ) 5*1 (~uj m ) (-w m ) 



-I * = *cq 



Inserting Eq. (46) we see that, due to the effective hop- 



ping amplitude , the second partial derivatives of the 
action still depend on the site distance i—j. This suggests 
to further apply a spatial Fourier transformation in order 
to simplify the calculations. Additionally, switching from 
Matsubara frequencies to continuous frequencies within 
a Wick rotation, yields the following explicit expression 
for Eq. (plj) 



= A(-oj, k) A* (u, k) - B* (-u, k) B(uj, k) 
with the abbreviations 



(52) 



A(w,k)=- 



1 



(0) 



(k, w ) 



- J(k) - - 



oi 0) (k;w,0|0,w)|* 



cq 



4 0) (k, w ) 



4 0) (k, ) 



1 2 ' 



and 
B(w,k)=- 



»i o) (k;0,0|-o;,w)|*eq| 



4 0) (k,o) 



,(n) 



(k, w )4 0) (k,- w ) 



(53) 



(54) 



Assuming a simple three-dimensional cubic lattice with 
lattice constant a yields J(k) = 2k Y^=i cos (^i a )- Tak- 



ing a closer look at Eq. ( 52 ) we see that it implicitly de- 



fines the dispersion relation uj (k). However, due to the 
complex expression found for af^ , a full evaluation of the 
above equations can only be done numerically. Just in 
the Mott insulator regime, where the 4-point correlation 
drops out due to <f cq = 0, it is possible to derive 
analytic expressions for the dispersions. In this special 

case we have to solve 1 = a^^k, uS) J(k) leading to the 
dispersion relations 



w±(k) 



± 



E, 



(«+!)- 



E, 



(n-l)- 



J(k) (tl 



(n+1)- 



{ 

(j(k) 



(n-1)- 



E(n+1)- - ^(k) 



t: 



±t 



(n+1) 



E, 



(n+l)- l n 



Ei 



E, 



E n - — E, 



(«+!)- 



)]} 



t 



(n+1) 



(n-1)- 

E n - 



E, 



(n+l)- 

(55) 



In order to clarify the physical meaning of equation ( 55 ) 
we set J(k) = 0, which yields the simple relations 



UJ h - 



— E„ 
= E, 



E, 



("+!)- 



E n - . 



(56) 
(57) 
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FIG. 2. (Color online) The pair excitation energy gap (left ) 
and the effective mass (right) in dependence of the temper- 
ature for the first lobe n = 1 at ^ e fl = Merit = 0.78<? with 
zero detuning A = and k — K cr it = 0.16g (solid), k = O.lg 
(dashed) and k = 10~ 4 g (dotted). The right picture repre- 
sents the effective mass of both the hole (lower branch/light) 
and particle (upper branch/dark) excitations. 



From these relations it is clear that Wh (w p ) is the 

energy needed to remove (add) a lower-branch polariton 
from a lattice site, which is occupied by n lower-branch 
polaritons. Therefore, we refer to these energies as the 
lower-branch-hole (h) and lower-branch-particle excita- 
tion (p), respectively. The same analysis can be applied 
to the full expression |52| ) leading to additional dispersion 
relation for upper-branch excitations Wh++ and w p ++ as 

well as mixed excitations Wh_| , LOh \-, w p _| and lj v |_. 

However, these additional excitation channels are hardly 
of any interest for our considerations because they occur 
at much higher energies compared to the lower-branch 
polariton excitations. The lowest dispersion relation in- 
volving upper polariton states is lj^ , which lies at an 

energy around 2g for n = 2. Therefore, we will just fo- 
cus in all further calculations on the lower-branch excita- 
tions as they completely determine the low-temperature 
physics of the system. In order to extract more detailed 
quantities of interest, we subsequently expand the dis- 
persion relations as follows 



w p ,h(k) 



ap 



2M 



p,h 



G(k 4 ). 



(58) 



From this expansion we are able to derive, within the 
Mott insulator phase, the energy gap E gap and the ef- 
fective mass of the particle and hole excitations M Pi h, 
respectively. As an example we evaluated their tempera- 
ture dependence for the first Mott lobe n = 1 with /Lt cr it 
and zero detuning in Fig. [2] We find that, both the pair 
excitation energy gap as well as the effective mass of par- 
ticle and hole excitations increase with higher tempera- 
tures. This effect becomes stronger as one approaches 
the critical hopping strength at the tip of the lobe. 
Moreover, we expect to find in the superfluid regime, 
apart from the gapped mode, also a linear excitation 
mode which is associated with the broken symmetry in 
the superfluid regime according to the Nambu-Goldstone 



theorem [55] . The dispersion relation for this mode reads 

w p , h (k)«c|k|+0(k 2 ). (59) 

The results obtained from a numerical evaluation of the 
above formulas is presented in Fig. [3] for a mean par- 
ticle density n = 2 at zero temperature and vanishing 
detuning A = 0. The pictures are arranged in a ta- 
ble in such a way that each line represents a specific 
physical quantity, for instance the second line Fig. |3(b)| 
shows plots for the energy gap. Furthermore, the plots 
in each column correspond to a fixed effective chemical 
potential. The left column shows plots for an effective 
chemical potential below the critical one, the middle col- 
umn represents the critical effective chemical potential, 
and the right column corresponds to an effective chem- 
ical potential above the critical one. In the first line in 
Fig. 3(a) we show the particle and hole excitation spec- 
tra in k = fc(l,l,l) direction. In the Mott phase and 
at the phase boundary we always observe two excitation 
modes corresponding to the particle (dotted) and hole 
(dashed) excitations, respectively. These modes are al- 
ways gapped in the Mott-insulator phase. By approach- 
ing the phase boundary at least one gap vanishes and 
thus for /x e g > /icrit the particle mode becomes gapless at 
the phase border, whereas for n B g < /j CI i t the hole mode 
becomes gapless at the phase border. When approaching 
the phase boundary exactly at the lobe tip both particle 
and hole modes become gapless. Going further into the 
superfluid regime we find a gapped excitation mode as 
well as the anticipated gapless linear mode. In Fig. |3(c) 
we plot the corresponding effective masses. In the Mott 
phase we observe the masses of both the particle (dotted) 
and hole (dashed) excitations, whereas in the superfluid 
phase also only one massive mode survives. Additionally, 
we depict the sound velocity of the polariton excitations 
in the superfluid phase in Figs. |3(d)| and |4j Figure |3(d)| 
shows the dependence of the sound velocity on the hop- 
ping strength. We find that it approaches a finite value at 
the tip of the Mott lobe, but vanishes at all other points of 
the Mott lobe. This behaviour shows that the JCH model 
has a dynamical critical exponent of z = 1 which has 
been recently confirmed in a large-scale quantum Monte- 
Carlo simulation by M. Hohenadler et al. [SB]. Entering 
the superfluid phase the sound velocity increases steadily. 
However, as pictured in Fig. [4] if the system is tuned out 
of resonance, i.e. A ^ 0, the sound velocity drops signif- 
icantly. Finally, we note that our results Figs. [3] and [4] 
are in good qualitative agreement with the results from 
Ref. [57]. However, due to the restriction to the lowest 
hopping order in the effective action, our results lose va- 
lidity deep in the superfluid phase. For this reason we 
do not obtain a shift of the maximum of the sound ve- 
locity in Fig. [4j as is observed in Ref. [57J ■ In order to 
obtain better results in the superfluid regime higher hop- 
ping corrections must be considered. Corresponding per- 
turbative approaches for higher order corrections for the 
JCH model have been numerically calculated in Ref. [55J . 
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SF modes: k z = 0.03 g 
MI modes: k z = 0.004 g 
Critical: kz = 0.0075 g 
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SF modes: kz = 0.015 g 
MI modes: kz = 0.01 g 
Critical: k z = 0.012 < 



15 




SF modes: k z = 0.02 g 
MI modes: kz = 0.005 g 
Critical: kz = 0.008 g 



(a) Particle (dotted) and hole (dashed) dispersion relations in the Mott phase (gray), on the phase boundary 
(black) and in the superfluid phase (solid) in the direction k = fc(l, 1, 1). 
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(b) Energy gap for particle (dotted) and hole (dashed) excitations in the Mott phase and for the massive 

mode in the superfluid phase (solid). 
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(c) Effective mass for particle (dotted) and hole (dashed) excitations in the Mott phase and for the massive 

mode in the superfluid phase (solid). 
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(d) Sound velocity for the massive mode in the superfluid phase. 



FIG. 3. (Color online) Various dynamic results for a mean particle density n = 2 at zero temperature and vanishing detuning 
A = 0. Left column at fj, e f{ — —OAg, middle column for tip of the lobe at jj, e g — ~0.37g and right column at fi e f{ = — 0.34gi. 
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FIG. 4. (Color online) Sound velocity c normalized by the 
sound velocity at zero detuning Co in dependence of the de- 
tuning parameter A for n = 2 with k = 0.00lK cr it (solid) and 
k = 0.15«: cr it (dashed). 



V. CONCLUSION 

In summary, we successfully applied the cumulant ex- 
pansion approach from Ref. |65j to derive a Ginzburg- 
Landau theory for the Jaynes-Cummings-Hubbard model 
up to fourth order in the symmetry-breaking currents and 
up to first order in the hopping strength. From the re- 
sulting effective action we extracted the phase diagram of 
the inherent quantum phase transition of the JCH model 
for finite temperature. In the case of vanishing temper- 
ature our result is in accordance with the results found 
from mean- field calculations [79l [82] ■ Subsequently, we 
derived the excitation spectra, energy gaps, and effective 
masses of the lower-branch polariton-excitations in the 
Mott insulator phase as well as in the superfluid phase. 
We investigated the temperature dependence of both the 
pair excitation energy gap and the effective mass of the 
particle and hole excitations in the Mott phase. Further- 
more, we analyzed how the sound velocity in the super- 
fluid phase depends at zero temperature on the hopping 
parameter and the detuning parameter. Finally, we point 
out that the the Ginzburg-Landau approach of this paper 
can be generalized to describe the real-time dynamics of 
the JCHM. This has already been shown for the real-time 
dynamics of the Bose-Hubbard model in Refs. [571 155] , 
We thank N. G. Berloff, C. Ciuti, and S. Schmidt for 
useful discussions and especially M. Hayn for valuable 
suggestions. 



Appendix A: Fourth Order Coefficient 



Here we evaluate the fourth order coefficient (33), 



which involves the expectation values of time ordered 
products of four operators. First we notice that, for the 
time-ordered product of two annihilation and two cre- 
ation operators, there are 6 distinct permutations leading 
to different expectation values. Each of these orderings 
itself has 4 time variable permutations corresponding to 



Ti <H- r 2 and r 3 ■<-» T4. Thus, overall one finds 24 terms 
for the expectation value. Luckily, the integrals over dif- 
ferent time-variable permutations yield the same result, 
and thus, they just lead to a fixed pre-factor 4. For this 
reason, one just needs to determine the 6 different ther- 
mal averages for one specific time-ordering. Furthermore, 
these expectation values are local quantities and, there- 
fore, we drop the site indexes in the following calcula- 
tions. Thus, one has to determine the following expres- 
sions: 

(a f (ri) cl>(t 3 ) a(r 2 ) a(r 4 )) , (a f (Vi) a(r 2 ) a\r 3 ) a(r 4 )) , 

(a(r4)a t (Ti)a t (T 3 )a(r 2 )) , (a(r 4 ) a(r 2 ) a f (Ti) a f (T 3 )) , 

(a(r 4 ) a^n) a(r 2 ) a f (T 3 )) , (^(n) a(r 4 ) a(r 2 ) a f (r 3 )) , 

(Al) 

With the help of the polariton mapping introduced in 
Ref. [7S], one can calculate these averages straightfor- 
wardly. Hence, we find for example the following expres- 
sion for the expectation value 



-. 00 



-/3E„ 



n=l 



P.7T=± 



x (B„ <s -B ( „_ 1 ,„)r 1 (E (n _ 1)7I -E np )r 2 , . 



x e 



(B„ p -B ( „-iv)r3 e {E( n -i )v -E na )Ti , , (AO) 



Subsequently, performing a Matsubara transformation 
according to ( 25 ) yields a formal integral of the form 



dte at / dt x e 
1) ^0 



1 = 1 
with the solution 
7=7 



Mi 



dt 2 e c 



dt 3 e dt3 (A3) 



(a+b+c+d)P 



1 



(a + b + c + d)(b + c + d)(c + d)d 

e (a+b+c)f3 _ 2 _ 1 _|_ p (a+b)l3 



(a + b + c)(b + c)cd b(a + b)c(c + d) 



-1 



ab(b + c)(b- 



d) 



(A4) 



The variables a, b, c, d correspond to differences of energy 
eigenvalues. Due to energy conservation these variables 
have to fulfill the condition a + b + c + d = and thus the 



above solution (A4) has a pole in the first term. There- 



fore, we have to take the limit 
lim 



e (a+b+c+d)/3 _ 1 

a+b+c+d^a (a + b + c + d)(b + c + d)(c + d)d 



> + c + d)(c + d)d' 



(A5) 



Since this pole arises for all expectation values (Al) we 



always have to consider this particular limit. Taking this 
result into account, the explicit expression for the expec- 



tation value ( A2 ) in Matsubara space is given by 
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1 



-l + /(-^l- E -l + "-A +g ">°) 

(^ ml +g_i + „ iA -£„, c ,)(-^ m i- J E_i + ^, A +£„, c ,)(^ m2 +g.i + „ iA -£ n , p ) 
Wm2 ^ W m 3 + £^_i +n ^ — E-i+ ntV 



+ 





(w m4 +-E-i+„ 


,„-B n ,c.)(-W m 4-B-l+r,,„+B„, c «)(-W m 3-B_i + „ i , 






W m 2 — W m 3 + E-i +7h \ — E-i+ n ^ v 










(w m2 +-E-i + „ 


.A-S„, p )(-cj ml +w m2 + S„, a -E„. p)(-W m3 -B-l + r 


,„+s„, P ) 



_ Wm3 + W m 4 — 



^n,a,A ^n,ot,v ^n,p,A ^n,p,v • 



(A6) 



However, this expression still possesses some poles for 
special choices of Matsubara frequencies. Fortunately, all 
these poles can be eliminated by investigating the corre- 
sponding limits analogous to equation ( A5 ) . Further care 



has to be taken considering the occurrence of the ground- 
state energy due to its uniqueness. Similar expressions 
can be calculated for the other expectation values ( Al I. 



A. Kay and J. K. Pachos, New J. Phys. 6, 126, (2004) 

P. Treutlein, T. Steinmetz, Y. Colombe, B. L. Lev, P. 

Hommelhoff et al., Fortschr. Phys. 54, 702 (2006) 

R. P. Feynman, Int. J. Th. Phys. 21, 467 (1982) 

H. S. J. van der Zant, F. C. Fritschy, W. J. Elion, L. 

J. Geerligs, and J. E. Mooij, Phys. Rev. Lett. 69, 2971 

(1992) 

A. van Oudenaarden and J. E. Mooij, Phys. Rev. Lett. 
76, 4947 (1996) 

M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. 

S. Fisher, Phys. Rev. B 40, 546 (1989) 

M. Naraschewski, H. Wallis, A. Schenzle, J. I. Cirac, and 

P. Zoller, Phys. Rev. A 54, 2185 (1996) 

M. R. Andrews, Science 275, 637 (1997) 

K. W. Madison, F. Chevy, W. Wohlleben, and J. Dal- 

ibard, Phys. Rev. Lett. 84, 806 (2000) 

J. R. Abo-Shaeer, C. Raman, J. M. Vogels, and W. Ket- 

terle, Science 292, 476 (2001) 

S. Kling and A. Pelster, Phys. Rev. A 76, 023609 (2007) 
S. Kling and A. Pelster, Laser Physics 19, 1072 (2009) 
D. M. Stamper-Kurn, M. R. Andrews, A. P. Chikkatur, 
S. Inouye, H.-J. Miesner, J. Stenger, and W. Ketterle, 
Phys. Rev. Lett. 80, 2027 (1998) 

F. Schreck, L. Khaykovich, K. Corwin, G. Ferrari, T. 
Bourdel et al., Phys. Rev. Lett. 87, 080403 (2001) 

S. Rothel and A. Pelster, Eur. Phys. J. B 59, 343 (2007) 
J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht et 
al. Nature 453, 891 (2008) 

G. Roati, C. D'Errico, L. Fallani, M. Fattori, C.Fort et 
al. Nature 453, 895 (2008) 

R. Graham and A. Pelster, Int. J. Bif. Chaos 19, 2745 
(2009) 

A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. 

Pfau, Phys. Rev. Lett. 94, 160401 (2005) 

A. R. P. Lima and A. Pelster, Phys. Rev. A 81, 021606 

(2010) 

T. van Zoest, N. Gaaloul, Y. Singh, H. Ahlers, W. Herr 
et al., Science 328, 5985 (2010) 



[22] D. Jaksch, C. Bruder, J.I. Cirac, C. W. Gardiner, and P. 

Zoller, Phys. Rev. Lett. 81, 3108 (1998) 
[23] M. Greiner, O. Mandel, T. Esslinger, T. W. Hansen, and 

I. Bloch, Nature 415, 6867 (2002) 
[24] D. Jaksch, H.-J. Briegel, J. I. Cirac, C. W. Gardiner, and 

P. Zoller, Phys. Rev. Lett. 82, 1975 (1999) 
[25] O. Mandel, M. Greiner, A. Widera, T. Rom, T. W. 

Hansen et al., Nature 425, 6961 (2008) 
[26] M. Riebe, H. Haffner, C. F. Roos, W. Hansel, J. Benhelm 

et al., Nature 429, 6993 (2004) 
[27] C. Roos, G. Lancaster, M. Riebe, H. Haffner, W. Hansel 

et al., Phys. Rev. Lett. 92, 220402 (2004) 
[28] L. Sanchez-Palencia and L. Santos, Phys. Rev. A 72, 

053607 (2005) 

[29] V. Ahufinger, L. Sanchez-Palencia, A. Kantian, A. San- 

pera, and M. Lewenstein, Phys. Rev. A 72, 063616 (2005) 
[30] D. Clement, A. Varon, M. Hugbart, J. Retter, P. Bouyer 

et al., Phys. Rev. Lett. 95, 170409 (2005) 
[31] K. V. Krutitsky, A. Pelster, and R. Graham, New J. 

Phys. 8, 187 (2006) 
[32] V. A. Yurovsky, A. Ben-Reuven, P.S. Julienne, and C.J. 

Williams, Phys. Rev. A 62, 043605 (2000) 
[33] S. J. J. M. F. Kokkelmans, H. M. J. Vissers, and B. J. 

Verhaar, Phys. Rev. A 63, 031601 (2001) 
[34] P. Wiirtz, T. Langen, T. Gericke, A. Koglbauer, and H. 

Ott, Phys. Rev. Lett. 103, 080404 (2009) 
[35] N. Gemelke, X. Zhang, C.-L. Hung, and C. Chin, Nature 

460, 995 (2009) 
[36] J. F. Sherson, C. Weitenberg, M. Endres, M.Cheneau, I. 

Bloch et al., Nature 467, 68 (2010) 
[37] W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon et al., 

Science 329, 547 (2010) 
[38] K. J. Vahala, Nature 424, 6950 (2003) 
[39] G. Khitrova, H. M. Gibbs, M. Kira, S. W. Koch, and A. 

Scherer, Nature Physics 2, 81 (2006) 
[40] T. M. Benson, S. V. Boriskina, P. Sewell, A. Vukovic, S. 

C. Greedy et al., Frontiers of Planar Lightwave Circuit 

Technology, pp. 39 (Springer, Berlin, 2005) 



13 



P. Grangier, G. Reymond, and N. Schlosser, Fortschr. [63 
Phys. 48, 859 (2000) 

M. Bayindir, B. Temelkuran, and E. Ozbay, Phys. Rev. [64 
B 61, 11855 (2000) 

M. J. Hartmann, F. G. S. L. Brandao, and M. B. Plenio, [65 
Nature Physics 2, 849 (2006) 

A. D. Greentree, C. Tahan, J. H. Cole, and L. C. L. [66 
Hollenberg, Nature Physics 2, 856 (2006) 

D. G. Angelakis, M. F. Santos, and S. Bose, Phys. Rev. [67 
A 76, 031805 (2007) 

M. Bajcsy, S. Hofferberth, V. Balic, T. Peyronel, M. [68 
Hafezi et al., Phys. Rev. Lett. 102, 203902 (2009) 
H. Deng, G. Weihs, C. Santori, J. Bloch, and Y. Ya- [69 
mamoto, Science 298, 199 (2002) 
J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. [70 
Jeambrun et al., Nature 443, 7110 (2006) 

R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West, [71 
Science 316, 1007 (2007) 

A. Amo, D. Sanvitto, F. P. Laussy, D. Ballarini, E. del [72 
Valle et al., Nature 457, 291 (2009) 
A. Amo, J. Lefrere, S. Pigeon, C. Adrados, C. Ciuti et [73 
al., Nature Physics 5, 805 (2009) 

A. Z. Muradyan and G. A. Muradyan, J. Phys. B 35, [74 
3995 (2002) 

J.-Q. Liao, Z. R. Gong, L. Zhou, Y.-X. Liu, C. P. Sun et [75 
al., Phys. Rev. A 81, 042304 (2010) 

T. Carmon, T. J. Kippenberg, L. Yang, H. Rokhsari, S. [76 
Spillane et al., Optics Express 13, 3558 (2005) 
F. P. Laussy, E. del Valle, A. Gonzalez- Tudela, E. Cancel- [77 
lieri, D. Sanvitto et al., 18th Int. Symp. "Nanostructures: 
Physics and Technology" 298 (2010) [78 
S. M. Spillane, T. J. Kippenberg, O. J. Painter, and K. [79 
J. Vahala, Phys. Rev. Lett. 91, 043902 (2003) [80; 
A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R. S. 
Huang et al., Nature 431, 162 (2004) 

T. J. Kippenberg, S. M. Spillane, D. K. Armani, and K. [81 
J. Vahala, Optics Letters 29, 1224 (2004) 

M. Trupke, J. Goldwin, B. Darquie, G. Dutier, S. Eriks- [82 
son, J. Ashmore, and E. A. Hinds, Phys. Rev. Lett. 99, 
63601 (2007) [83 
M. Bajcsy, S. Hofferberth, V. Balic, T. Peyronel, M. 
Hafezi et al., Phys. Rev. Lett. 102, 203902 (2009) [84 
F. S. F. Brossard, X. L. Xu, D. A. Williams, M. Hadji- 
panayi, M. Hopkinson et al., Appl. Phys. Lett. 97, 111101 [85 
(2010) [86 
A. Laucht, J. M. Villas-Boas, S. Stobbe, N. Hauke, F. 
Hofbauer et al., Phys. Rev. B 82, 075305 (2010) [87 

[88; 



G. Ctistis, A. Hartsuiker, E. van der Pol, J. Claudon, W. 
L. Vos et al., Phys. Rev. B 82, 195330 (2010) 

F. Nissen, S. Schmidt, M. Biondi, G. Blatter, H. E. 
Tiireci, and J. Keeling, arXiv: 1202. 1961 

B. Bradlyn and F. E. A. dos Santos, Phys. Rev. A 79, 
013615 (2009) 

F. E. A. dos Santos and A. Pelster, Phys. Rev. A 79, 
013614 (2009) 

T. D. Grass, F. E. A. dos Santos, and A. Pelster, Laser 
Physics 21, 1459 (2011) 

T. D. Grass, F. E. A. dos Santos, and A. Pelster, Phys. 
Rev. A 84, 013613 (2011) 

A. Hoffmann and A. Pelster, Phys. Rev. A 79, 053623 
(2009) 

D. McKay and B. DeMarco, Rep. Prog. Phys. 74, 054401 
(2011) 

D. F. Walls and G. J. Milburn, Quantum Optics 
(Springer, Berlin, 1995) 

E. T. Jaynes and F. W. Cummings, P. IEEE 51, 89 
(1963) 

S. Trotzky, L. Pollet, F. Gerbier, U. Schnorrberger, I. 
Bloch et al., Nature Physics 6, 998 (2010) 

F. Gerbier, S. Trotzky, S. Foelling, U. Schnorrberger, J. 
D. Thompson et al., Phys. Rev. Lett. 101, 155303 (2008) 
S. Sachdcv, Quantum Phase Transitions (Cambridge 
University Press, Cambridge, 1999) 

L. D. Landau and E. M. Lifschitz, Course of Theoretical 

Physics, Volume 5 (Akademie Verlag, Berlin, 1970) 

J. W. Negele and H. Orland, Quantum Many-Particle 

Systems (Westview Press, Boulder, 1998) 

W. Metzner, Phys. Rev. B 43, 8549 (1991) 

J. Koch and K. Le Hur, Phys. Rev. A 80, 023811 (2009) 

H. Kleinert and V. Schulte-Frohlinde, Critical Properties 
of 4> 4 Theories (World Scientific Publishing Company, 
Singapore, 2001) 

J. Zinn-Justin, Quantum Field Theory and Critical Phe- 
nomena (Oxford University Press, Oxford, 2002) 
S. Schmidt and G. Blatter, Phys. Rev. Lett. 103, 086403 
(2009) 

P. Buonsante and A. Vezzani, Phys. Rev. A 70, 033608 
(2004) 

I. Carusotto, D. Gerace, H. E. Tureci, S. De Liberato, C. 
Ciuti et al., Phys. Rev. Lett. 103, 033601 (2009) 

J. Goldstone, Phys. Rev. 127, 965 (1962) 

M. Hohenadler, M. Aichhorn, S. Schmidt, and L. Pollet, 

Phys. Rev. A 84, 041608 (2011) 

S. Schmidt and G. Blatter, Phys. Rev. Lett. 104, 216402 
(2010) 

C. Hei and W. von der Linden, arXiv: 1105 . 2418 



